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Abstract 

In multibody simulation, the Gear-Gupta-Leimkuhler method for only persistent contacts enforces constraints 
on position and velocity level at the same time. It yields a robust numerical discretization of differential alge¬ 
braic equations avoiding the drift-off effect. In this work, we carry over these benefits to impacting mechanical 
systems with unilateral constraints. For this kind of a mechanical system, adding the position level constraint 
to a timestepping scheme on velocity level even maintains physical consistency of the impulsive discretiza¬ 
tion. Hence, we propose a timestepping scheme based on Moreau’s midpoint rule which enables to achieve not 
only compliance of the impact law but also of the non-penetration constraint. The choice of a decoupled and 
consecutive evaluation of the respective constraints can be interpreted as a not energy-consistent projection to 
the non-penetration constraint at the end of each time step. It is the implicit coupling of position and velocity 
level which yields satisfactory results. An implicit evaluation of the right hand side improves stability prop¬ 
erties without additional cost. With the prox function formulation, the overall set of nonsmooth equations is 
solved by a Newton scheme. Results from simulations of a slider-crank mechanism with unilateral constraints 
demonstrate the capability of our approach. 
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1 Introduction 


Dynamical motion with impacts plays an important role in the characterization of general mechanical systems 
at least after discretization in space. The monographs |171|2l[TTl|T4l[T7l summarize the state-of-the-art physical, 
mathematical and numerical setting of this kind of impacting mechanical systems: 
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if gy < 0, then 0 < gt -f egj ± Ay > 0 . 
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Starting from the initial conditions ([Hi-©, the development of the system’s state given by position q and veloc¬ 
ity V is described by a non-impulsive behavior (HJl almost everywhere. It is influenced by the generalized mass 
matrix M and right hand side forces h. Due to the Signorini-Moreau condition closed or opening sclero- 
nomic contact gaps g affect this type of motion by varying contact force parameters A. The notation g -L A 
stands for g^A = 0. The force parameters weight the columns of in the equations of motion (@1). For count¬ 
able time instances tj, the velocities jump enforced by an impact Ay according to ©. Newton’s impact law d?]) 
provides the respective relationship of active pre- and post-impact velocities using the kinematic coefficient of 
restitution e. With dfl]), the position q can be calculated by the fundamental theorem of calculus for weakly 
differentiable functions. 



Figure 1: Slider-crank mechanism with unilateral constraints. 

One might think that the differential complementarity problem ([Tll-(17]l is integrated best by an event-driven 
time-integration strategy. However as event-driven schemes resolve the exact time of impact, they cannot 
consistently model Zeno phenomena, i.e. infinite impacts occurring in a finite time interval. Timestepping 
schemes discretize the equations of motion including the constraints consistently without resolving the exact 
transition points. Robustness benefits from this approach but the accuracy is comparably low. 

1.1 Moreau’s midpoint rule 

Because of the consistent approach within timestepping methods,, we focus on a well-established representative, 
i.e. Moreau’s midpoint rule lIT^ . It summarizes both impulsive and non-impulsive phases: first by calculating 
-in some sense- the mean impulsive force within fixed time steps At and second by incorporating the results 








implicitly in a time-discretization on velocity level: 
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Equation (fTOl) approximates the positions by their explicitly calculated midpoint values and the velocities by 
the respective explicit values at the beginning of the time step. The only unknown variables in (I8]l-(l9ll are the 
positions and velocities v„+i at the end of the time step and the discrete mean impulsive force A„+i. With 
the explicit predictor 


8 m = S {qn + y 


< 0 , 


( 11 ) 


an active set of constraints indicated by the subscript ( )red is determined and every corresponding interaction 
can be treated by Newton’s impact law on velocity level. Next to classical impacts also closed contacts, opening 
contacts or impacts occurring without a normal contact impulse are naturally included. The complementarity 
formulation of Newton’s impact law can be equivalently reformulated by dint of the prox function related to 
a convex set C C IR; this is often easier to solve numerically than a complementarity problem ifTSl . As the 
proximal point proX(-(x) € C of x G IR is defined as 


proX(-(x) = argmin ||x —x*|| , (12) 

.v*GC 

Newton’s discrete impact law corresponds row-by-row to 

An+l,red = PrOXj]^+ (A„_|_i j-gd “ ^ (kn+l,red “h ^k«,red)) • (1^) 

This expression can be solved iteratively with the -in the easiest case- positive diagonal parameter matrix r 
controlling the speed of convergence; in this work, we choose a Newton scheme without adapting r as solution 
method. As termination criterion, the natural monotony test or a tolerance for the residuum can be employed. 


1.2 Problem statement 

Moreau’s midpoint rule holds for general impacting mechanical systems. We reveal improvement possibilities 
at a glance of a planar impacting slider-crank mechanism H. 


1.2.1 Slider-crank mechanism with unilateral constraints 


For the slider-crank mechanism in Figure [T] angles q = { 61 , 62 , 63 )^ and angular velocities v = (ft)i,ft) 2 ,G) 3 )^ 
rely on an absolute description concerning an inertial x-y-frame of reference. The crank (1) has mass mi, 
rotational inertia around the center of gravity Ji and length /i. The connecting rod (2) is similarly represented 
by 1712, Ji and h- The slider (3) with m 3 and J 3 has height 2b and length 2a. Its center of gravity (x 3 ,y 3 ) is not 
fixed on one y-posifion buf can move wifhin a nofch of heighf d and clearance c. The gap functions are defined 
as illusfrafed in Figured 
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Figure 2: Definition of the gap functions for the slider-crank mechanism with unilateral constraints. 


The tangential gap functions are neglected because the frictionless case is considered. Accordingly, the con¬ 
straint matrix satisfies 

( —/i cos 01 —ZicosOi /i cos 01 Zicos0i \ 

—Z2COS02 —Z2COS02 Z 2 COS 02 Z 2 COS 02 | . (18) 

acos 03 -|-Zjsin 03 —a cos 03 -|-Zj sin 03 —a cos 03 -|-Z? sin 03 a cos 03 -|-Z? sin 03 / 

Assuming gravifation g in negative y-direction, the example fits exactly in the concept of a general impacting 
mechanical system ([I])-©: 

■/i (x + "^2 + «I3) ZiZ2COs(0i -02) (^-hm3) 0\ 

ZiZ2Cos(0i-02)(f+m3) + + 0 , (19) 

0 0 73/ 

-ZiZ2sin(0i -02) (^-hm3)a)|-gZiCOS0i +m2 + m3)\ 

ZiZ2sin(0i-02) (^-hm3)a)f-gZ2Cos02(^-hm3) j . (20) 

1.2.2 Simulation results 

With the characteristics of Bl reprinted in Table [H we analyze the movement of the center of gravity of the 
slider (3) using Moreau’s midpoint rule and the time step size At = 10^^ s. Figure [3] shows the results for 
four different coefficients of restitution being the same for each contact possibility, respectively. The curves 
presented here, as well as the graphs depicted in the original literature ||4l, show the violation of the non¬ 
penetration condition ®, especially in the simulations with a low coefficient of restitution. In Figure [3] (a), the 
center of gravity of the slider (3) exceeds 10^^ m and in the later course it falls below the value —10^^ m, what 
corresponds to a pervasion of the slider (3) with the bordering wall. Also for higher coefficients of restitution, 
the gap functions fall below zero, but of course for shorter time periods. 

Figure |4] shows the time curve of the gap functions and their time derivatives for a coefficient of restitution 
£ = 0.1. The gap between the slider (3) and the surrounding wall is small as well as the initial configuration 
03 „ = 0 and the support at the center of gravity prevent the revolution of the slider (3). Hence, pairs of gap 
functions on opposite sides appear almost symmetric. Obviously, the non-penetration condition is violated. A 
detailed view shows the drift-off effect: the gap drifts approximately linearly to negative values while the gap 
velocity is slightly negative. However, drift does not have a dominant effect for this configuration because the 
negative gap functions remain comparatively small in contrast to the geometric dimensions. When contacts stay 
closed for longer time periods, the drift-off effect will not be negligible anymore. 
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1.2.3 Outline 


For systems with only persistent contacts, one could consider the constraints on position level. For systems 
which undergo impacts in addition, this would yield non-consistent discretizations. Hence, when both impacts 
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Figure 3: Movement of the center of gravity of the slider (3) for different coefficients of restitution. 


and longer periods of closed contacts occur, neither of these two approaches, i.e. neither on position nor veloc¬ 
ity level seems to be satisfactory. Gear, Gupta and Leimkuhler @ proposed a solution to a related problem for 
only persistent contacts: they enforce constraints on position and velocity level at the same time. The additional 
constraint equation is compensated by a second set of Lagrange multipliers. The purpose of this work is to apply 
the Gear-Gupta-Leimkuhler method to systems with unilateral constraints to overcome the drift-off effect for 
closed contacts as roughly indicated in fll. Thereby, we summarize and extent our student work ifT^ . First, we 
present the Gear-Gupta-Leimkuhler method for a slider-crank mechanism without clearance. The application 
to unilateral contacts shows that a decoupled strategy satisfying velocity and position level constraints one after 
the other is not energy-consistent. A unified formulafion which fakes info accounf fhe impacf law as well as 
fhe non-penefrafion consfrainf af fhe same lime lurns oul lo be a successful approach. We close fhe paper wifh 
some open questions for fulure research direcfions. 
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Figure 4: Gap functions and their time derivatives for £ = 0.1. 


Geometrical characteristics 

h = 0.1530m 
h = 0.3060m 
a = 0.0500m 
b = 0.0250m 

c = 0 . 0010 m 

Inertia properties 

nil = 0.0380 kg 
m2 = 0.0380 kg 


m 3 = 0.0760kg 

Ji = 7.4- lO^^kgm^ 

72 = 5.9-10-^kgm2 

73 = 2.7- lO-^kgm^ 

Force elements 

g = 9.81m/s^ 

Initial conditions 

010 = 0.0 

02 o = 0-0 

03o = 0-0 
coio = 150.01/s 
ft) 2 „ = -75.01/s 
= 0 . 01 /s 


Table 1: Characteristics of the slider-crank mechanism with unilateral constraints. 























































2 Gear-Gupta-Leimkuhler method for persistent contacts 


We explain the expected effect of the Gear-Gupta-Leimkuhler method |l 6 l on the numerical solution of a unilat¬ 
erally constrained mechanical system with the help of a bilaterally constrained slider-crank mechanism adapted 
from im. The drift-off effect is analyzed for constraint formulations on position, velocity and acceleration level 
as well as for the Gear-Gupta-Leimkuhler formulation Q. 

2.1 Slider-crank mechanism with bilateral constraints 

Figure |5] shows a bilaterally constrained slider-crank mechanism. The angles q= ( 61 , 62 )^ are chosen as gen- 



Figure 5: Slider-crank mechanism with bilateral constraints. 


eralized coordinates, the angular velocities v = (cOijCOz)^ as generalized velocities. With the same notation as 
in Section [U we gain equations of motion which are more specific than stated in ([Tll-(17]l: unilateral contacts 
condense to bilateral constraints and impacts never occur. The generalized mass matrix satisfies 


M = 


_f Ji-' 1 -)n2 Z1Z2COS (01 — 62) +WI3) 

ZiZ2COs(ei-02)(f+^3) /2 + /|(f+'”3) 


and fhe vector of generalized forces is given by 

^ /-ZiZ 2 sin( 0 i - 02 ) {^ + mi,) (O^-gh cos di +m 2 + m 2 ) 
Zi Z 2 sin (01 - 02 ) (^ -h m 3 ) mf - gZ 2 cos 02 (^ -h m 3 ) 

The bilaferal consfrainf holds fhe slider (3) af a fixed y-posifion 

g = h sin 01 -h h sin 02 = 0 

by causing a consfrainf force in direction of 




Zi cos 01 
I 2 cos 02 


( 21 ) 


( 22 ) 


(23) 


(24) 


2.2 Simulation results 

The simulafions are accomplished wifh fhe characferisfics of Table[T]and fhe time sfep size At = 10^'^ s. A direcf 
compufafion of fhe consfrainf compliance considering fhe consfrainf g on posifion level yields a differenlial 
algebraic system of index 3. If is known to be badly condifioned and e.g. scaling of fhe consfrainf equation 
yields an heuristic improvemenf concerning fhe sfabilify of fhe numerical infegrafion scheme ||9||. Instead of 
thaf, we focus on replacing fhe consfrainf by ifs respecfive fime derivatives, which improves fhe robusfness of 
numerical solvers consisfenlly from an analytic poinf of view. Arnold 131 mentions fhis sfrafegy in fhe confexf 
of index reduction. A draw back of index reduction is fhe driff-off effecf 0. Figure 0 shows fhe roughly 
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Figure 6: Drift-off effect of the slider (3) for equations formulated on velocity and acceleration level. 


linear development of the y-position of the slider (3) for a long-time simulation of the index 2 system, i.e. 
considering the constraint g on velocity level. Figure [6] also displays the drift-off effect for the constraint 
formulation on acceleration level g, i.e. for the index 1 system. As presented in Q, the drift-off is expected 
to be parabolic and in fact the y-position increases with rising gradient. As a compromise of both robust 
simulation and asymptotically little drift-off, one usually considers the constraints on velocity level. To even 
overcome the linear drift-off effect in the index 2 system. Gear, Gupta and Leimkuhler proposed a formulation 
which considers the constraints on position and velocity level simultaneously m. The original index 2 system 
extends to 


q = v+W'^Y, (25) 

Mv = h + W^)i, (26) 

g = 0, (27) 

g = 0. (28) 


The Lagrange multiplier if/ compensates the added equation and the constraint is satisfied on position as well 
as on velocity level maintaining the stability of the index 2 formulation. 








3 Gear-Gupta-Leimkuhler method for unilateral contacts 


We analyze two extensions of Moreau’s midpoint rule (cf. Section [TTI) . A decoupled approach turns out not to 
be energy-consistent. A unified approach meets our expectations but demands the computational effort of an 
implicit solution scheme. 


3.1 Decoupled approach 


The adaption of Moreau’s midpoint rule is performed by adding a correction term enforcing the non-penetration 
constraint at the end of each time step: 


Vn+\ —Vn+MfJ -|-W^A„+i) , (29) 

A,)+l,red = PJ'®*IR+ (An+yred “ ^ jej “h ^Snjed)) i (^0) 

^n+1 = “I 2 , (31) 

'¥„+i = proXiR+ ('T„+i - rg„^i) . (32) 


As in Moreau’s midpoint rule, the calculation of the velocities v„+i is achieved by using the average Lagrange 
multiplier A„+i. This computation is decoupled from the calculation of the positions with the average 
Lagrange multiplier Consecutively, the velocities v„+i are used to determine an explicit forecast: 


9n+l 


9n + 


Vn+l+Vn 
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At . 


(33) 


It is used as an initial value for the iterative computation of and 

The positive effect of low computational effort, due to a decoupled calculation of the two different vectors 
of Lagrange multipliers, is subtended by the low physical accuracy of the results. This is clarified by fhe 
developmenf of fhe entire energy confenf of fhe slider-crank mechanism wifh unilateral constraints using the 
characteristics of Table [T] (cf. Figure [TOl) . No energy sources are applied but for a coefficient of restitution 
e =0.1, our simulation results with time step size At = 10^^ s reveal a fluctuating entire energy content and do 
not show the expected decreasing trend. Hence, the timestepping scheme (|2^- (|3^ does not provide a valid and 
physical accurate model of a system underlying unilateral constraints. 


3.2 Unified approach 

What is the problem in (I29l) - (l32l) ? The equations of motion are derived via an energy principle, the impact law 
results from Newton’s admittedly kinematic considerations. However, the additional term does not 

correspond to any physical principle but can be interpreted as part of the Karush-Kuhn-Tucker conditions for a 
projection at the end of each time step concerning the Euclidean metric: 

min||^„+i, (34) 

g(^„+i)>0. (35) 

It is not astonishing that the entire energy content oscillates. Studer ifTSl mentions this type of discretization 
discussing the bouncing ball example. As a workaround, it is suggested to introduce a penetration tolerance 
depending on the specific setting. To our opinion, fhe ferm has fo be coupled wifh (|2^ fo ensure a 

physical accurafe behavior in general. Our proposition is presenfed in fhe nexf section. 


3.2.1 Discretization scheme 

A possible coupling is the implicit evaluation of the constraint matrix Wm = ^ ^ ^ maintaining the nice 

properties of the midpoint concept, e.g. symplecticity iH. As this strategy already enforces the solution of a 





_ L ( gn+l+gn Vn+l+Vr. 


nonlinear system of equations, we additionally evaluate the generalized force vector Hm = h 
implicitly to benefit from a more stable discretization of important stiffness conttibutions. The generalized mass 
matrix comprises geomettic nonlinearities and its implicit evaluation needs comparatively large effort. Hence, 
an explicit evaluation Mm is chosen: 
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2 ’ 2 

Adding the active constraints and defining a nonlinear system of equations, only the dependency on the un 
known variables 


is interesting: 
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In contrast to Section [H the discretizations Hm = /im (^«+i)V«+i) and Wm = ^M(^n+i) explicitly depend on 
the unknown values and v„+i. The roots of the reduced system of equations (xred) can be solved using 
Newton’s method 


m+1 — yin _ 

■^«,red -^n.red 


dx. 
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With the scleronomic gap functions being simplified concerning effective evaluations 

§«+l =Wn+iqn+l -g„ + WM (qn+l) At+W M (^,,+ 1) (^n+l) i 
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the derivative of iPred (-^red) with respect to Xred can be deduced by eliminating rows and columns corresponding 
to inactive contacts from the following matrix 
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With an appropriate function /, a distinction of cases is necessary for the rows containing the prox function: 
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ProxiR+ (/(x 
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Analyzing exemplary the more difficult if-case of the prox function, the derivatives of the third row of (p are 
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Specify geometry, start time t = 0, end time tE, time step size At, parameters and initial configuration 


For t < tE 

Evaluate Mm (flOl) 


Evaluate gM and define active sef (fTTI) 


Evaluafe g and g 


While Newfon iferafion nof converged 

Compute Ji’m, W'm, g"*, g'”, proX]R+1^,„ and (p‘ 


Compute ^ 
Update (|4^ 


1 dWM 1 



Ix^’ dx lx™’ 

dx lx"” 

Sx 1 y" 




and 


3(p I 
dx Ijc"' 


§3 


Wrife resulf of fime step 


Update t by At 


Write result 


Eigure 7: Elowchart of the proposed unified timestepping scheme. 


Eor the fourth row, it is 
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As long as the active set is empty, all Eagrange multipliers are equal to zero and the system of equations is 
solved without the constraint part. As soon as the active set is not empty, the system of equations is extended 
by Newton’s impact law and the non-penetration constraint, respectively. The algorithm is set up as shown in 
Eigure |7J 


3.2.2 Simulation results 

In Eigure m the results concerning the slider-crank mechanism with unilateral constraints and characteristics as 
in Tableware presented using a time step size At = 10^^ s. The qualitative behavior is similar to the behavior for 
Moreau’s midpoint rule shown in Eigure[3l Especially for high coefficients of restitution, the patterns resemble. 
The change in the theoretical framework mainly affects persistent contacts, which rarely occur for e > 0.5. 
In contrast for e =0.1, the drift-off effect has a comparatively high influence. The proposed timestepping 
scheme yields a distinct change in the system’s behavior. The drift-off effect does no longer occur and the 
non-penetration condition is satisfied improving the physical accuracy in comparison to Eigure [3l 

The development of the gap functions and their time derivatives for £ = 0.1 is presented in Eigure |9l The 
drift-off effect has vanished and the gap functions are not negative anymore. The gap velocities are still slightly 
smaller than zero in time periods where the drift-off effect occurred in the previous simulations. However, this 
slow trend to increasing permeation is compensated by the second set of Eagrange multipliers enforcing the 
non-penetration constraint. 

To further investigate the physical accuracy of the method, the qualitative development of the entire energy 
content for Moreau’s midpoint rule and the unified Gear-Gupta-Eeimkuhler approach is shown in Eigure [TOl 
The entire energy content after four seconds of simulation differs slightly. As Moreau’s midpoint rule does 
not ensure the compliance of the constraints, a reference line is shown based on the same algorithm as the 
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Figure 8: Movement of the center of gravity of the slider (3) for different coefficients of restitution. 


unified Gear-Gupfa-Leimkuhler approach buf only enforcing fhe impacf law like Moreau’s midpoinf rule and 
neglecfing fhe non-penefrafion condition. The reference line is based on (l3^ wifhouf fhe lasf row and wifhouf 
fhe ferm (^n+i)'F„+i in fhe firsf row. The energy developmenf for fhe reference system is slightly smaller 
than for the unified Gear-Gupfa-Leimkuhler approach and a bif higher fhan for Moreau’s midpoinf rule. The 
proposed approach leads fo fhe same qualifafive behavior of fhe entire energy confenf and fo slighfly differenl 
quanfifative resulfs. 

Due fo fhe implicif discretization, fhe computing time increases by a factor of fen in confrasf fo fhe explicif 
Moreau’s midpoinf rule when using fhe same lime step size. However, we gain a sfable discretization which 
allows comparatively larger fime step sizes for sfiff problem formulalions. 
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Figure 9: Gap functions and their time derivatives for £ = 0.1. 
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Figure 10: Energy content for different formulations. 























































4 Conclusion 


Within this work, we propose a timestepping scheme for impacting mechanical systems with unilateral con¬ 
straints. The new scheme is based on Moreau’s midpoint rule and enables to achieve not only compliance of the 
impact law but also of the non-penetration constraint. It is shown that the decoupled application of the Gear- 
Gupta-Leimkuhler method for bilateral constraints can be interpreted as a projection to the non-penetration 
constraint at the end of each time step. As this strategy does not lead to an energy-consistent discretization, our 
proposition couples position and velocity level with an implicit evaluation of the constraint matrix in the frame¬ 
work of midpoint discretizations. Without significant additional cost, we also approximate the right hand side 
in the same manner to achieve enhanced stability properties. Adding the active constraints in each time-step 
by means of the prox function concept leads to a system of nonsmooth equations which is solved by a Newton 
scheme. 

Results from simulations of a slider-crank mechanism with unilateral constraints demonstrate the overcoming of 
the drift-off effect and the performance of our unified approach. If is reduced concerning the differential index 
and insofar better conditioned than position level discretizations. Concerning impacting mechanical systems, it 
is even more important that our scheme is physically consistent due to the impulsive concept and the implicitly 
incorporated projection. Because of the implicit discretization, the computation time increases significantly in 
comparison to Moreau’s midpoint rule using the same non-controlled time step size. However, for stiff problem 
formulations, our proposition should result in a more stable discretization and possible larger time step size 
choices. 

We do not have applied our scheme to a stiff problem. An analysis concerning numerical issues could be ad¬ 
dressed by utilizing backward error analysis ||3|. Using this concept, the interpretation of the induced projection 
on the non-penetration constraints should be discussed in addition. The relationship to the kinetic metric and 
appropriate extensions should be studied ifT^ . Finally, the consideration of friction taking into account the 
overdetermined differential algebraic setting according to ifTOl would extend our work. 
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